function b = adsecurity_gr(Xit,Pspot,MU,CAPH,T,rho)
%This function gives the ad security for the future period given history.
%The output is b=[b(lambda=1|Lambda_t); b(lambda=2|Lambda_t),...,b(lambda=n|Lambda_t)]

%In the first step we get the current period, spot price and guarantted
%price from history Lambda_t
    t0=length(Xit)+1;
    mu0=MU(:,:,t0);
    tildeP=findtildeP(Xit,Pspot,MU);
    pspott=Pspot(:,t0);
    n=length(mu0);
    tildeP=min([pspott repmat(tildeP,n)],[],2);
    %MU = repmat(MU_ins_t,[1 1 T]);
    if t0==T
        %if current period is final period, 
        % The bond is spot price minus min(spot price, guarnateed price)
        b=mu0-tildeP;
        b(b<0)=0;
    else
        %If not, we first get all future spot price.
        tf=T-t0;
        Pf=Pspot(:,(t0+1):T);
        % future costs
        muf = MU(:,:,t0+1:T);

    
        %Then we calcualte the aggregate probability agent would not lapse
        %(same as denominator of equation 1 in the paper)
        %First, we get the state where the agent would stay
    
        Nstates = size(CAPH,1);
        Pf=reshape(Pf,[Nstates 1 tf]);
        stay  = repmat(permute(Pf,[2,1,3]),Nstates,1)>=repmat(tildeP,1,Nstates,tf);
    
        % death is a no stay state
        stay(:,Nstates,:) = 0;
        % probability of getting to each state in a 
        % given period -- transit from periods where i stay 
        capH=CAPH(:,:,t0);
        pt = zeros(Nstates,Nstates,tf);
        pt(:,:,1) = capH;
        
        E = zeros(Nstates,Nstates,tf);
    
        E(:,:,1) = rho*pt(:,:,1).*stay(:,:,1).*repmat(permute(muf(:,:,1),[2,1,3]),Nstates,1,1);

        D = zeros(Nstates,Nstates,tf);
        D(:,:,1)= rho*pt(:,:,1).*stay(:,:,1);
        for t=2:tf
            capH=CAPH(:,:,t0+t-1);
            pt(:,:,t) = pt(:,:,t-1).*stay(:,:,t-1)*capH;
            E(:,:,t) = rho^t*pt(:,:,t).*stay(:,:,t).*repmat(permute(muf(:,:,t),[2,1,3]),Nstates,1,1);
            D(:,:,t)= rho^t*pt(:,:,t).*stay(:,:,t);
        end
        % Add them up

        % numerator
        E = sum(sum(E,2),3);
        num = mu0+E;
        %denumerator
        den = 1 + sum(sum(D,2),3);
        %use equation 3 in ra note to calculate bond.

        b=num-den.*tildeP;
        %b(b<0)=0;
    end
end